import numpy as np


def f(x):
    '''
    定义f(x)
    :param x: 需要计算的函数值
    :return: 返回sinx/x 的值
    '''
    if x == 0:
        x = 1e-20
    return 1 / x


def romberg(a, b, e):
    '''
    使用龙贝格算法，计算sinx/x在(a,b)之间的差距
    :param a: 积分下限
    :param b: 积分上限
    :param e: 误差极限
    :return:
    '''
    h = b - a
    # n为划分的区间数
    n = 1
    # k表示行数
    k = 0
    # T为龙贝格积分表,定义一个20*20 的矩阵。矩阵的值为0
    T = np.zeros((15, 15))
    # T_0=就等于梯形面积
    T[0][0] = (f(a) + f(b)) * (h / 2)
    # 错误率
    err = 10
    result = 0
    while err >= e:
        # S表示梯形的面积
        S = 0
        # 行数+1
        k += 1
        # h为[a,b]这段区间中，被分成n分，每一份的长度。
        # 假设a=5,b=7,n为5分,h为[5,7]这段区间，被分成了5份，每份0.4,h=0.4
        # h = (b - a) /n   h=h/2,因为后面n=n*2
        h = h / 2
        # 变步长梯形法求第0列 Tn,n越大
        for i in range(1, n + 1):
            S = S + f(a + (2 * i - 1) * h)

        # T_2n=1/2T_n+h/2*积分和
        # T[k][0]放置的是T_{2k}的值，第0列,T_0,T_2,T_4,T_8
        T[k][0] = T[k - 1][0] / 2 + h * S
        # 计算S_n,C_n,R_n
        for m in range(1, k + 1):
            T[k][m] = T[k][m - 1] + (T[k][m - 1] - T[k - 1][m - 1]) / (4 ** m - 1)

        # 为了增加梯形的精度，所以增加划分区间的
        n = n * 2
        # 计算误差率
        err = abs(T[k][k] - T[k - 1][k - 1])
        result = T[k][k]

    # np.savetxt("result.txt", T, delimiter=',', fmt="%.10f")
    return result


# 初始化初始值，a,b,e
# a表示积分下限
# b表示积分上限
# e表示误差限
a = 1
b = 3
e = 1e-5
print("龙贝格算法答案计算1-3的积分值---", romberg(a, b, e))
